rm(list = ls())  #清理环境
#setwd("C:/Users/mortal/Documents/ASQ")  #设置工作路径

#加载所需要的包
# install.packages(readxl)  # 如已安装，可注释该步
library(readxl)    # 读取速度优于 read.xlsx()，且不依赖 Java
library(xlsx)      # 输出数据时要用到 write.xlsx()
library(geepack)   # gee分析
library(mice)      # 缺失值处理
library(stringr)   # 处理字符串

#数据读取
data <- read_excel("./output_data1.xlsx",sheet = 3) #读取第二个工作簿“数据1”
str(data)                         #查看数据结构
colnames(data)                    #查看变量名
#data_com=data[complete.cases(data[,1:5]),]把这句加入循环，去获得非空变量

#孕周转小数
a <- as.numeric(str_sub(data$week,1,2))
b <- as.numeric(str_sub(data$week,4,4))
a[is.na(a)]<-0
b[is.na(b)]<-0
data$week <- a+b/7
data$week[data$week==0]<-NA

#孕周重编码
data$week_f[data$week<37 |data$week>42]<-2
data$week_f[data$week>=37 & data$week<=42]<-1

data_test <-data      #原本删除了不分析的变量，19和34
length(unique(data_test$id))      #查看样本量
############################处理缺失值##########################
#转化数据类型
vars1 <- colnames(data_test)[c(1:43)]
data_test[vars1]<- lapply(data_test[vars1],as.numeric)

#查看数据的缺失情况
#vars2 <- colnames(data_test)[c(1,14:49)]
#index1 <- is.na(data_test[vars2])
#rowSums(index1==T)  #查看研究对象数据缺失情况（每行观测值的缺失数量）
#colSums(index1==T)  #查看变量缺失情况（每列变量的缺失数量）
#data1 <- data_test[rowSums(index1==T)>15,]，剔除超过15个变量缺失的观察对象

### 转化数据类型
colnames(data_test) 
vars3 <- colnames(data_test)[c(2:8,26,29,30,33,35,39,40)]
vars4 <- colnames(data_test)[c(9:25,27,28,31:32,34,36:38,41:45)]
data_test[vars3]<- lapply(data_test[vars3],as.numeric)
data_test[vars4]<- lapply(data_test[vars4],factor)
str(data_test)

'###缺失值填补
#截取暴露变量并去重
vars2 <- colnames(data_test)[c(1,14:49)]
data_mice <- data_test[vars2]
data_mice <- unique(data_mice)
str(data_mice)

data_mice$id <- as.character(data_mice$id) 
#这一步不能少，mice默认填补方法是“pmm”,ID如果是数值变量回导致数据矩阵奇异

miss<- md.pattern(data_mice)   #是显示缺失的数据是什么模式的
attach(data_mice)
imp <- mice(data_mice,m=5,seed = 1)#m是表示要生成的完整数据集的数量method是默认pmm
fit <- with(imp,glm(diversity~fruit+egg+meat+beans+potato,
                    family = binomial(link = "logit")))
#理论上这里with()没必要，比较with（）当中不是gee模型，可直接complete()
pooled <- pool(fit)
data_mice1 <- complete(imp,action = which.min(pooled$glanced$AIC))  
#选择 AIC 与BIC最小模型的填补数据集？意义
detach(data_mice)

write.xlsx(data_mice1,"data_mice1.xlsx")
#用于拼接的ID、月份与ASQ结果，merge全连接自动扩充
data_asq<-data_test[,c(1:13)]
data_mice1$id <- as.numeric(data_mice1$id)
data_comp<-merge(data_asq,data_mice1)

write.xlsx(data_comp,"data_comp.xlsx")'



#######################ASQ结局（二分类）与各自变量的单因素分析##############
data_anal <- data_test
#data_anal <- data_anal[,-1]
### 转化数据类型
'vars3 <- colnames(data_anal)[c(1:13,33:34,37,39,43,46)]
vars4 <- colnames(data_anal)[c(14:32,35:36,38,40:42,44:45,47:49)]
data_anal[vars3]<- lapply(data_anal[vars3],as.numeric)
data_anal[vars4]<- lapply(data_anal[vars4],factor)'
str(data_anal)
colnames(data_anal)

#拆分自变量与因变量
varsx <- colnames(data_anal)[c(13:27)]
varsy <- colnames(data_anal)[c(9:13)]
#varsz<-colnames(data_anal)[c(8:12)]#先单因素吧，协变量没整理好
data_x <- data_anal[varsx]
data_y <- data_anal[varsy]
str(data_y)
#data_y=as.list(data_y)

###for循环将分析结果输出
for (j in 1:length(colnames(data_y))) {
  data_run=data_anal[complete.cases(data_x[,1]),]
  varsx <- colnames(data_run)[c(14:25,27)]
  varsy <- colnames(data_run)[c(9:13)]
  data_xrun <- data_run[varsx]
  data_yrun <- data_run[varsy]
  output<- data_yrun[,j]
  input<- data_xrun[,1]
  fit <- geeglm(formula = output~input,
                data = data_anal,
                id=id,family = "binomial",
                corstr = "ar1", scale.fix=TRUE)
  a <- summary(fit)
  out<-a$coefficients
  out$variable <- colnames(data_x)[1]
  for (i in 2:length(colnames(data_x))) {
    data_run=data_anal(complete.cases(data_x[,i]),)
    varsx <- colnames(data_run)[c(14:25,27)]
    varsy <- colnames(data_run)[c(9:13)]
    data_xrun <- data_run[varsx]
    data_yrun <- data_run[varsy]
    input<- data_xrun[,i]
    output<- data_yrun[,j]
    fit <- geeglm(formula = output~input,
                  data = data_anal,
                  id=id,family = "binomial",
                  corstr = "ar1", scale.fix=TRUE)
    a <- summary(fit)
    out1<-a$coefficients
    out1$variable <- colnames(data_x)[i]
    out <- rbind(out,out1)  #rowbind加行
  }
  write.xlsx(x=out,file = "outcom2.xlsx",sheetName = colnames(data_y)[j],append = T)
}

###修改时间 2021/9/15
